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ABSTRACT 

We develop a new formalism to study the dynamics of fluid polytropes in three dimen- 
sions. The stars are modeled as compressible ellipsoids and the hydrodynamic equations 
are reduced to a set of ordinary differential equations for the evolution of the principal axes 
and other global quantities. Both viscous dissipation and the gravitational radiation reac- 
tion are incorporated. We establish the validity of our approximations and demonstrate 
the simplicity and power of the method by rederiving a number of known results concerning 
the stability and dynamical oscillations of rapidly rotating polytropes. In particular, we 
present a generalization to compressible fluids of Chandrasekhar's classical results for the 
secular and dynamical instabilities of incompressible Maclaurin spheroids. We also present 
several applications of our method to astrophysical problems of great current interest such 
as the tidal disruption of a star by a massive black hole, the coalescence of compact binaries 
driven by the emission of gravitational waves, and the development of instabilities in close 
binary systems. 
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1. INTRODUCTION 



In a series of recent papers (Lai, Rasio & Shapiro 1993a,b, 1994a,b, hereafter LRS1-4 
or collectively LRS), we have developed a new analytic method to calculate ellipsoidal 
figures of equilibrium both for a single rotating polytrope and for polytropes in binary 
systems. Our approach is based on the use of an energy variational principle to construct 
approximate equilibrium solutions and to determine their dynamical and secular stability. 
In addition to providing compressible generalizations for all the classical incompressible 
solutions (Chandrasekhar 1969, hereafter Ch69), our method can also be applied to more 
general models of binary systems where the stellar masses, radii, spins, entropies, and 
polytropic indices are all allowed to vary independently for each component. As a result, 
we were able to study the equilibrium and stability properties of many different types of 
fairly realistic binary models (see especially LRS4). 

We have used the method to explore the implications of instabilities for the coalescence 
of close binary systems (LRS2, LRS3). Of particular importance currently are coalescing 
neutron star binaries, which are the primary targets for the detection of gravitational 
waves by LIGO (Abramovici et al. 1992; LRS3). In our previous papers we have studied 
these systems by considering equilibrium sequences of binary configurations, treating the 
binary separation r as a dynamical variable but with all other variables assuming their 
equilibrium values. This is adequate to capture the essential features of the evolution, but 
does not provide quantitatively accurate results when the orbital evolution takes place on a 
timescale comparable to the internal hydrodynamic timescale. The purpose of the present 
paper is to extend our formalism and develop a fully dynamical theory for the evolution 
of compressible ellipsoids. 

Detailed studies of hydrodynamic interactions between stars normally require large- 
scale numerical simulations in three dimensions (see, e.g., Lai, Rasio, & Shapiro 1993c; 
Rasio & Shapiro 1991, 1994). These simulations demand extensive computational resources 
and cannot be used to cover a large parameter space. In our treatment, we replace the 
infinite number of degrees of freedom in a fluid system by a small number of variables 
specifying the essential geometric and kinematic properties of the system. The dynamics 
is then described approximately by a set of ordinary differential equations (ODEs) for the 
time evolution of these variables. This represents an enormous simplification. However, 
for many problems, useful insights and even reliable quantitative results can be obtained 
using such an approach. The simplicity of a formulation in terms of ODEs allows for 
extensive coverage of the parameter space. In addition, we can follow the evolution of a 
system over a large number of dynamical times without having to worry about excessive 
computational time or about the growth of numerical errors. This advantage is crucial for 
studying the secular evolution of a system on a dissipative timescale while still allowing 
for dynamical processes. A distinctive example is the coalescence of binary neutron stars, 
which begins as a very slow orbital decay driven by the emission of gravitational waves, 
but ends in a rapid hydrodynamic merging of the two stars after the stability limit has 
been reached (LRS3; Rasio & Shapiro 1994). 

In the incompressible limit, the dynamics of an isolated ellipsoid was first formulated 
by Lebovitz (1966; see also Ch69, Chap. 4). This so-called Riemann-Lebovitz formulation 
of the problem forms the basis of many subsequent studies and astrophysical applications 
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of ellipsoidal models. In particular, the effects of fluid viscosity and gravitational radiation 
reaction on the stability of a rapidly rotating star were investigated by Chandrasekhar 
(1970), Press & Teukolsky (1973), Miller (1974), and Detweiler & Lindblom (1977). The 
Riemann-Lebovitz equations have also been used to study the gravitational collapse of 
rotating stars and the resulting gravitational radiation by relaxing the assumption of in- 
compressibility and incorporating a polytropic equation of state while still assuming a 
homogeneous density profile (see Shapiro 1979 and references therein). Nduka (1971) has 
studied the dynamics of an incompressible ellipsoid in a parabolic orbit around a fixed 
point mass. 

For compressible systems, Carter & Luminet (1983, 1986) have developed a dynam- 
ical theory of an affine stellar model in the context of tidal interactions with a massive 
black hole. Kochanek (1992a,b) has used the method to study the tidal capture of a poly- 
trope by a point mass and the coalescence of two polytropes in a close binary. In the 
affine model, the fluid compressibility is treated under the assumption that the surfaces 
of constant density remain self-similar ellipsoids. Our dynamical model for compressible 
ellipsoids is essentially equivalent to the affine model (cf. LRS1), but our formulation of the 
problem is quite different and our dynamical equations more closely resemble the Riemann- 
Lebovitz equations in the incompressible limit. In contrast to previous studies, we present 
a completely general formulation for isolated stars and for a star on a bound or unbound 
trajectory around a point mass. We incorporate both viscosity and gravitational radiation 
reaction as possible dissipation mechanisms. Our formulation makes explicit use of global 
quantities such as the total angular momentum and fluid circulation which are conserved 
in the absence of dissipation. This greatly simplifies the description of many dynamical 
processes, and we feel that this simplification has not been fully appreciated in previous 
studies based on the affine model. 

The main purpose of this paper is to formulate general dynamical equations for com- 
pressible ellipsoids and to present a small survey of astrophysical applications. We wish to 
illustrate how these general equations can be used to tackle a variety of multidimensional 
hydrodynamic problems of great current interest. The dynamical equations for an isolated 
ellipsoid supported by a polytropic equation of state are derived in §2, where our general 
assumptions are also presented (§2.1). Many of the expressions derived in §2 can be read- 
ily transported to more general situations. In §3 we consider the dynamical stability and 
oscillations of a single rotating polytrope. In §4, we incorporate viscous dissipation and 
gravitational radiation reaction into the dynamical equations, and we study the secular 
evolution of isolated rotating stars. In §5, the dynamical equations derived in §2 for single 
stars are generalized to a binary system (either bound or unbound) consisting of a poly- 
trope and a point mass. The dynamical stability of general Roche-Riemann binary models 
is considered in §6. In §7, we study tidal encounters of a star with a massive body, and we 
compare our results with those of linear perturbation theory as well as previous nonlinear 
calculations. In §8, we consider the evolution of binary systems driven by viscous dissipa- 
tion. In §9, we incorporate gravitational radiation reaction in our treatment of binaries, 
and we study the orbital evolution driven by gravitational radiation. 
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2. DYNAMICAL EQUATIONS FOR SINGLE STARS: 
RIEMANN-S ELLIPSOIDS 



In this section, we derive the equations of motion for an isolated compressible Riemann- 
S ellipsoid. Our equations represent a generalization of the Riemann-Lebovitz equations, 
which only apply to an incompressible fluid (Ch69, Chap. 4). 

2.1 Assumptions 

Throughout this paper, we adopt a polytropic equation of state for the fluid, 

P = Kp 1+1 / n . (2.1) 

Our basic assumptions concerning the interior structure of the star can be summarized as 
follows (see LRS1 for more details). We assume that the surfaces of constant density can 
be represented approximately by self-similar ellipsoids. The geometry is then completely 
specified by the three principle axes of the outer surface. Furthermore, we assume that the 
density profile p(m) inside each star, where m is the mass interior to an isodensity surface, 
is identical to that of a spherical polytrope with the same volume. In the reference frame 
comoving with the star's center of mass, the velocity field of the fluid is modeled as a 
linear superposition of three components: (1) a rigid rotation of the ellipsoidal figure; (2) 
an internal fluid circulation with uniform vorticity; and (3) an ellipsoidal expansion or 
contraction (see eqs. [2.3] and [2.4] below). 

Under these assumptions, the number of internal degrees of freedom for each star is 
reduced to nine in general: the three principal axes a\, a^-, 03 of the outer surface, the 
three components of the angular velocity of the ellipsoidal figure fii, O27 ^3, and the 
three components of the vorticity (j, £2, C3- We shall further restrict ourselves to the 
cases where both the angular velocity and the vorticity are parallel to one of the principal 
axes, here taken to be the z or 03 axis by convention. This choice corresponds to the 
Riemann ellipsoids of type-S (Ch69). As we show in §2.2, under these conditions, the 
exact hydrodynamic equations (the Euler and Navier-Stokes equations, with and without 
gravitational radiation reaction) are replaced by a set of ODEs for the time evolution of 
the principle axes ai, a2, 03, the angular velocity of the ellipsoidal figure O = Q3, and the 
angular frequency of the internal circulation A oc (3 (cf. LRS1, §5.1). 

As we noted in the introduction, our model is formally equivalent to the affine model 
developed extensively by Carter and Luminet (1985; Luminet & Carter 1986). In the 
incompressible limit (n = 0), both models can be obtained equivalently by imposing fixed 
holonomic constraints on the system, requiring the fluid surface to be ellipsoidal and the 
fluid velocity to be a linear function of coordinates. This is known as the Dirichlet problem 
(see Ch69, Chap. 4). In the n = limit, our dynamical equations reduce explicitly to the 
equations obtained in the Riemann-Lebovitz formulation of the Dirichlet problem (Ch69, 
§27). For a single isolated star with n = 0, the solutions we derive represent exact solutions 
of the true hydrodynamic equations. For n^O, our solutions are only approximate since 
the isodensity surfaces can no longer be exactly ellipsoidal, and the velocity field of the 
fluid cannot be described exactly by a linear function of coordinates (see Ipser & Managan 
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1981 for a formal proof). For binary systems, our solutions are always approximate because 
we truncate the gravitational interaction to quadrupole order (see §5 below). 

Clearly, our formalism cannot be applied to hydrodynamic processes such as stellar 
collisions that are violent enough to disrupt a star completely or to produce significant 
shock heating. Even for dynamical processes that are not very disruptive, our treatment 
allows us to follow only a small subset of the the many degrees of freedom that may be 
important. This is especially true for compressible systems, which have typically many 
more important degrees of freedom than incompressible systems. Consider for example 
small nonradial oscillations of a star (see, e.g., Cox 1980). For n = 0, only the f-modes 
of oscillation exist. But when n ^ 0, there are many additional p-modes and possibly 
also g-modes of oscillation, corresponding to nonellipsoidal disturbances of the star. Such 
nonellipsoidal motions cannot be described at all within our simplified treatment. Even 
for n = 0, an ellipsoidal model can only represent the / = 2, quadrupolar f-modes of 
oscillation. These modes often dominate the response of a star to tidal perturbations 
by passing objects (see Kochanek 1992a and §7 below), but they are not sufficient for a 
complete description of the dynamical response (cf. Press & Teukolsky 1977). 

2.2 Derivation of The Dynamical Equations 

Let ei, e2, and e3 be the basis unit vectors along the instantaneous directions of the 
principal axes of the ellipsoid, with e3 along the rotation axis (we refer to this as the "body 
frame"). In the inertial frame, the velocity field u inside the ellipsoid can be written as 

u = u s + u e . (2.2) 

Here u s is the velocity field of an equilibrium Riemann-S ellipsoid, 

u s = (^-A-Qjx 2 e 1 + (-^A + Qjx^, (2.3) 

where O is the angular velocity of the ellipsoidal figure, and A is the angular frequency of 
the internal fluid circulation (cf. LRS1, §5.1). For nonequilibrium configurations, we add 
the velocity u e describing an ellipsoidal expansion or contraction of the star, 

ai a 2 a 3 

u e = — xiei H x 2 e 2 H x 3 e 3 . (2.4) 

a\ a 2 a 3 

The kinetic energy is then given by 

T = Jd 3 x^pu 2 =T S +T e . (2.5) 

Here T s is the "spin" kinetic energy (rotation and internal circulation of the fluid; cf. LRS1, 
eq. [5.6]), 

T s = -/(A 2 + O 2 ) - -K n M ai a 2 An, (2.6) 
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where I is the moment of inertia, 

1= ^« n M(a? + a|), (2.7) 
o 

K n is a structure constant depending on the polytropic index (numerical values are given in 
LRS1, Table 1), and T e is the kinetic energy associated with the expansion or contraction, 

T e = ^K n M(al + a 2 2 + al). (2.8) 

The total internal energy U of the fluid is given by 

U = J^-dm = k 1 Kp 1 J n M 1 (2.9) 

where k\ is a constant depending on the polytropic index, and p c is the central density, 
equal to that of a spherical polytrope with the same mass and volume in our approximation 
(see LRS1, §2.1). The self-gravitational potential energy is given by 

w = - JL°£/, (2 , ) 

5 — n R 

where R = (aia 2 a 3 ) 1 / s is the mean radius of the ellipsoid, and 

/=2§2> with I = A ia j + A 2 a 2 2 + A 3 a 2 3 (2.11) 

(/ = 1 for spherical configurations). The dimensionless index symbols A t are defined as in 
Ch69 (§17). 

The Lagrangian governing the dynamics of the ellipsoid can now be written as 

L( qi ;qi) = T-U-W, (2.12) 

where {qi} = {ai,a2,a3,(f),ip}, and {g^} = {di, a 2 , a 3 , A} (we have introduced angles 
(j) and ifj such that = and ip = A). The dynamical equations are obtained from the 
Euler-Lagrange equations 

-— = — . (2.13) 
dt d<ji dqi 

Using the relation (Ch69, Chap. 2) 



equation (2.13) for qi = a± can be written as 

-K n Md 1 = -K n Mai(0 2 + A 2 ) - -ft n Ma 2 OA + b^LEl <2 ^—MGpa 1 A u (2.15) 

5 5 5 nai p c 5 — n 
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where p = 3M/ (4ttR 3 ) is the mean density. Equations for a 2 and a 3 can be similarly 
obtained. For q t = and tp, equation (2.13) yields the conservation laws for angular 
momentum and circulation, 



dJ, 



dt 



0, 



dC 
~dl 



0, 



where J s is the angular momentum, 



dL 2 
J s = tttt = Ift - -K n Maia 2 A, 
oil 5 

and C is the fluid circulation (cf. LRS1) 

dL 2 
C = — = IA - -K n Maia 2 tt. 

OA 5 

The complete set of dynamical equations can be rewritten in the form 



ai = a i (0 2 + A / ) -2a 2 OA 



a 2 = a 2 (fl +A 2 ) -2aiOA 
2ttG 



a 3 



-a 3 A 3 p + 



5/Ci \ P c 1 



ra«„ y p c ai 
5/ci \P C 1 



nK n ; p c a 2 



«Kn / Pc a 3 



— (aiO — a 2 A) = — diO + d 2 A, 

(J/L 



d_ 

dt 



(— a 2 + aiA) = d 2 fi — diA, 



(2.16) 



(2.17) 



(2.18) 



(2.19) 

(2.20) 

(2.21) 

(2.22) 
(2.23) 



where q n = K n (l — n/5). Rather than using equations (2.22) and (2.23) as such, we shall 
use the equivalent pair 



A = 



a 2 


ai 


ai 


a 2 


a 2 


ai 


ai 


a 2 



O A 



a 2 ai 

O A 
— + — 

ai a 2 



2 — + — di-2 — + — d 2 



O A 



di - 2 



ai a 2 

O A 
— + — 

a 2 ai 



a 2 



(2.24) 
(2.25) 



In addition to the conservation of total angular momentum and circulation (eqs. [2.17]- 
[2.18]), it is easy to check that the total energy is also conserved, i.e., 



£ = S s =T + U + W = constant. 



(2.26) 



For di = 0, one can verify that equations (2.19)-(2.21) reduce to the equilibrium equations 
for Riemann-S ellipsoids derived in LRS1 (§5.1). 



7 



2.3 The Pressure Term 

To implement the dynamical equations in numerical integrations, it is convenient 
to express the pressure term (5kiP c ) / (n,K n p c ) in terms of _R , M and other dynamical 
variables, where R Q is the radius of a nonrotating spherical polytrope of the same mass M 
and polytropic constant K. This is done as follows. Since P c / p c = Kp X J n oc R~ 3 / n , we 
can write 

^l^ = C n ^-3/n ; (227) 
TW^ra Pc 

where C n is a constant depending only on n, K, M. We can obtain the value of C n by con- 
sidering a single, nonrotating spherical polytrope in equilibrium. For such a configuration, 
equations (2.19)-(2.21) become 







2tvG 



Qn tio 

As A\ = 2/3 for a sphere, we get 



R A lp - + C n R- 3 / n ^-. (2.28) 



Therefore, in equations (2.19)-(2.21), we can use 

5/ci P c GM (R ^ 3/1 



C n = —&J n -\ (2.29) 



nn n p c q n R Q \ R 



(n ^ 0). (2.30) 



Expression (2.30) is obviously not valid in the incompressible case. For n — > 0, ki/n — > 
2/5 and K n — > 1, so that 

M±Pc = 2Pc (n = Q) (231) 

nn n p c p c 

To evaluate 2P C / p c in this limit, we need to use the incompressible condition 

aia2(i3 = Rl = constant, (2.32) 



which gives 



^ + ^ + ^ = 0. (2.33) 
ai 02 as 



Now adding (2.19)/ai, (2.20)/a2 and (2.21)/a3, and using the above expression, we get 
5fci P c 2P C 



nn n p c p c 

-l 



(n = 0) 
(2.34) 



where we have used A\ + A 2 + A 3 = 2. From equation (2.30), we see that the internal 
energy (2.9) can also be written as 

R n GM 2 //0 3/n 



U = k 1 —M=- — . (2.35) 

p c 5 — n R Q \ R 



Note that U = for n = since h = (cf. LRS1, eq. [2.9]). 



3. DYNAMICAL OSCILLATIONS AND STABILITY OF SINGLE STARS 

In this section, we use the dynamical equations derived in §2 to study the stability 
and dynamical oscillations of Riemann-S ellipsoids. We first establish a general criterion 
for linear stability to small-amplitude dynamical perturbations (§3.1). We then study 
the dynamical oscillations of an ellipsoid, including finite-amplitude self-similar pulsations 
(§3.2) and small-amplitude nonradial oscillations (§3.3). Most of the results obtained in 
this section are well-established. The derivations presented here will demonstrate that our 
formalism can reproduce these results or generalize them to the compressible case. 



3.1 Dynamical Stability Criterion 

Here we show explicitly that the onset of dynamical instability along an equilibrium 
sequence can be determined from an appropriate energy function defined along that se- 
quence, as we have done in our previous papers (LRS) . To do so we only need to recast the 
dynamical equations into Hamiltonian form. From the Lagrangian (2.12), the canonical 
momenta associated with the generalized coordinates {qi} = {ai,a 2 ,a3,(f),ip} are calcu- 
lated as Pi = dL/dqi, giving 

P a% = -K n Ma h P^ = J s , P^ = C. (3.1) 
5 



The Hamiltonian is then 



H 



2(k u M/5) 



(p z +p 



1 +P 2 )+E(a u a 27 a 3 ;J sl C)- 



(3.2) 



The energy function E plays the role of an effective potential, and is given by (cf. LRS4, 
eqs. [2.23] and [2.24], where we make explicit the conserved quantities J s and C) 

E( ai , a 2 , a 3 ; J s , C) = ^-{J s + Cf + ^-(J s - Cf + U + W, (3.3) 

where 

I ± = l Kn M( ai Ta 2 ) 2 . (3.4) 
o 
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Hamilton's equations can then be written as 

1 dE 

-K n M'di = - — , (3.5) 
5 ddi 

together with J s = C = 0. For linear perturbations we write 

a>i = a ijeq (l + cti), (3.6) 
where |aij| <C 1, and equation (3.5) becomes 

d 2 E 



dciidcij 



didj 



ctj (no sum over i), (3.7) 



3 - - ■" 

where Uj = (K n /5)Maf 5ij. If we let aij oc e J<Tt , equation (3.7) gives 

Here a subscript "eq" means evaluating the function for equilibrium configurations. Equa- 
tion (3.8) is the condition that determines the onset of instability along a sequence of 
equilibrium configurations (see LRS1, §2). 

3.2 Homologous Pulsations 

Consider the homologous perturbations described by 

°>i(t) = f (*H, e <i, (3-9) 

about a compressible Riemann-S equilibrium configuration. Here £(£) is not necessarily 
close to unity (i.e., we allow for finite-amplitude oscillations). Such oscillations only exist 
for a compressible (n 7^ 0) configuration. Our description of the oscillations by equa- 
tion (3.9) is only approximate (cf. §3.3), but it is a good approximation for slowly rotating 
configurations. 

Adding (K n Ma 1 /5)x eq. (2.19), (/j n Mo 2 /5) x eq. (2.20), and (K n Ma 3 /5) x eq. (2.21) 
we obtain 

It ^ = 2T ,J2M + (?hMf,) r a/. (3 . 10 ) 
C \ np c J eq 

where I t = In + 122 + -^33, T s is given by equation (2.6), and we have used equation (2.10) 
and the fact that P c j p c oc (a^as) -1 ^ 71 . Since T s can be written as (LRS4, eqs. [2.23]- 
[2.24]) 

Ts = ^(J s +CY + ±-(J s -C)\ (3.11) 
and J s , C remain constant during a dynamical perturbation, we have 

T s = T s , eq C 2 - (3.12) 
10 



For the last term in equation (3.10), we use equation (2.30) and the expression 



R 

R 



eq 



f 1-2- 



W 



-n/ (3— n) 



(3.13) 



(LRS1, eq. [4.25]). Equation (3.10) then reduces finally to 

iti = iwi(r 1_3/n - r 2 ) + 2T s (r 3 - r 1 " 37 "), 



(3.14) 



where we have suppressed the subscript "eq". This equation was first derived for an 
axisymmetric, uniformly rotating polytrope by Tassoul (1970). We see that it applies also 
to our more general compressible Riemann-S ellipsoids. 

Consider now small-amplitude oscillations, with £ = l+£i and |£i| <C 1. Linearization 
of equation (3.14) yields 



\w\ 



2T 

(4-sr)--^(5-3r) 



(3.15) 



We see that stability requires 



2T S /\W\ 



3 S(1-2T S /\W\Y 
The corresponding oscillation frequency is given by 



(3.16) 



\W\ 



(3r-4)-2^(3r-5) 



(3.17) 



This is a well-known result (the Ledoux formula), first derived for uniformly and slowly 
rotating stars (see Tassoul 1978, §14.2). We see that it can also be applied (approximately) 
to more general ellipsoidal configurations, even when the rotation rate is large. 

For nonrotating stars, equation (3.17) gives the fundamental pulsation mode frequency 



\W\ 



(3r - 4) 



GM 1 

Rl Qn 



(3r-4), 



(3.18) 



where q n = (l—n/5)K n , and K n is defined by equation (2.7). The values of K n for different n 
are given in Table 1 of LRS1. In Table 1 here, we compare the predictions of equation (3.18) 
with the exact results obtained by integrating numerically the radial pulsation equation 
for polytropes with V = Ti (Cox 1980, Chap. 8). We see that equation (3.18) predicts the 
exact results remarkably well. The discrepancy remains less than ~ 3% for all values of 
n. For n = 3 and in the limit as n — * 0, equation (3.18) yields the exact results. This 
is not surprising since, in these limits, small radial pulsations in the fundamental mode 
are indeed homologous. The excellent agreement for all n gives us confidence that our 
treatment of small dynamical perturbations is always a very good approximation. 
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We can also consider the stability to finite-amplitude oscillations (Tassoul 1970). In- 
tegrating equation (3.14) once, we obtain 



hti 2 + \W\*(£) = constant, (3.19) 

where 

*m = T,l\W\-j 1 - 2TJ\W\ 

K > e 3(r-i){ 3 ( r -»' ( 1 

From equation (3.19), we see that \I/(£) acts like an effective potential for the oscillation. 
It is easy to see that £ = 1 is a local minimum of \I/(£) only when the condition (3.16) is 
satisfied. But even then it is not necessarily a global minimum. Consider the behavior of 
*(£) for £ » 1, 



1-2T S /|W/| 
3(r-l)£ 3r " 4 



(3.21) 



We see that when T < 4/3, \&(£) has a global minimum at £ = oo. Thus for T crit < T < 4/3, 
the star is only metastable (Tassoul 1970). The condition for stability to finite-amplitude 
dynamical perturbations is V > 4/3, independent of rotation. 

3.3 General Oscillations of a Maclaurin Spheroid 

We now consider the general linear oscillations 

ai(t) = a ijeq (l + ati), (3.22) 

with |a!j| <C 1. For simplicity we assume cii jeq = a2, e? , i-e., the equilibrium configuration 
is a compressible Maclaurin spheroid. 

We must linearize equations (2.19)-(2.23) to obtain the dynamical equations for c^. In 
the derivation, one must be careful to note that a Maclaurin spheroid need not be assigned 
the angular velocity of the frame in which it is at rest, i.e., the value of Q eq and A eq must 
be determined, like the frequency cr, from the analysis (Rossner 1967). Here we will use 
equation (3.7) directly. For a Maclaurin spheroid we have J s = —C (cf. eqs. [2.17]-[2.18] 
for cii = 02), and the kinetic energy term is simply 

J2 

Ts = ( Kn /5)M(ai + a 2 ) 2 ' (3 ' 23) 
Setting the first derivatives of E equal to zero, we obtain the equilibrium conditions 

T s = -E(a?^ - a 2 3 A 3 ), 

1 2 „ ( 3 - 24 ) 

n 



where we have defined 



3GM 2 

W^W' (3 ' 26) 
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so that W = EX. For convenience, we also define 



£i >= { aia >^k) eq - (3 - 26) 
After some algebra we obtain the following expressions, 

£11 =£ 22 = -(l + -)u- a 2 S(3B n - 2A X ) + \t 8 , 
n \ n J I 

£33 = - (l + - V + 2a 3^13, 

n V n I 



\lJ -ajZ(B 11 -A 1 ) + ^ 



(3.27) 













Us / 



£13 = £23 = \u- a 2 E(B 13 - A x ), 

where the Bij are defined as in Ch69 (Chap. 3) and the right-hand sides are evaluated for 
the equilibrium configuration (the subscript "eq" has been dropped). Let ai oc e irjt . The 
linearized equations for small oscillations then give 

£11 — luc 2 £12 £13 

£12 £11 - /11a 2 £13 J ( « 2 ] = 0. (3.28) 

£l3 £l3 £33 — -^330" 2 

We now determine the eigenfrequencies and eigenmodes by solving the linear sys- 
tem (3.28). One of the eigenvalues can be obtained from 

I lia 2 = Sn - S 12 = -U- a 2 E(2B n - A x ). (3.29) 

n 

Using the equilibrium conditions (3.24), we obtain 

a 2 = — Bii-O 2 , (3.30) 

Qn 

where a = a/inGp) 1 / 2 , O = 0/(yrGp) 1/2 , and p = 3M/(4irR s ) is the mean density. The 
corresponding eigenmode has the form 



(3.31) 



This mode is the compressible generalization of the toroidal mode found in Ch69 (or bar 
mode). In a reference frame corotating with the unperturbed star, fluid elements oscillate 
at a frequency a Q given by 

/ 4 \ V2 

a = tt±a = tt±l —Bn - O 2 . (3.32) 

\Qn J 

13 



fa x \ 






«2 


M 




\«3 I 







This expression agrees with the result (5.50) of Ch69 for n = 0. We see that the onset of 
dynamical instability is given by 



& = & dyn = -fl u . (3.33) 
The toroidal mode is neutrally stable as seen in the corotating frame (a Q = 0) when 



2 



Cl 2 = n 2 bif = —Bn, (3.34) 



corresponding to the point where the Jacobi sequence bifurcates from the Maclaurin se- 
quence (Ch69; LRS1). 

The eigenfrequencies for the two other modes are determined from 

(£n + S12 ~ hicr 2 )(S 33 - I 33 a 2 ) = 2£ 2 3 , (3.35) 

and the corresponding eigenmodes have the form 



(3.36) 









h i 


\a 3 J 


\a 3 /ai / 



For these two modes the linearized equations reduce to 

Sn + S 12 - hia 2 £ 13 \ ( a x 



2£l3 ^33 - ^33^ / V «3 



(3.37) 



For n^O, these two modes are somewhat similar to the homologous pulsations discussed in 
§3.2. They are also the compressible generalization of the zonal modes discussed by Ch69 
and Tassoul (1978). Since a 3 ^ oti, the deformation is not homologous and the results 
obtained in §3.2 are slightly modified. In particular, setting a = in equation (3.37), we 
obtain the critical F for dynamical stability, 



, = 4 2T S 
crit 3 3IW1 



1- 2 + 



A, \ T. 



W 13 -2A 3 J \W 



- s 



(3.38) 



This stability condition is slightly more restrictive than equation (3.16) since a more general 
trial function has been used here. With some algebra it can be shown that expression (3.38) 
agrees with the result derived in LRS1 (eq. [6.3]). 

The expressions for the eigenfrequencies are quite complicated for general n. For n = 
the result is simpler. Subtracting the two equations in (3.37) from each other, taking the 
n — > limit, and using the equilibrium conditions (3.24), we obtain 



2 1 Wai 

Ilia 2 + (4a?Bn - 2a§B 13 )E (3aiS 13 - 2a|A 3 )E - J 33 a 2 U 3 



(3.39) 
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a 2 = 



i+ff, (» = 0). (3.40) 



The eigenfrequency is given by 

4B 11 + ^§(6B 33 -4B 13 
This is identical to the result given by Ch69 (eq. [5.63]). The corresponding eigenmode is 

(ra = 0). (3.41) 



«1 ) 








H 




«3 y 




V-2 



For a nonrotating (spherical) star with n = 0, the toroidal mode (3.31) is simply the / = 2, 
m = ±2 Kelvin mode, while the mode (3.41) describes axisymmetric pulsations with / = 2, 
m = 0. With By = 4/15 for a sphere, equations (3.30) and (3.40) both give a 2 = 16/15, 
an exact result for the frequency of the Kelvin mode. 

Note that the "transverse-shear modes" discussed by Ch69 (§33) cannot be incorpo- 
rated in our treatment, since they involve perturbations of the rotation axis. 

4. DISSIPATIVE EFFECTS FOR A SINGLE STAR 

We now incorporate dissipative forces into the dynamical equations derived in §2. In 
general, dissipation modifies the Euler-Lagrange equations according to 

d dL dL ^ , , „ . 

dt dqi dqi 

where T qi is the generalized force associated with the coordinate qi. The generalized 
dissipative forces are defined so that the dissipation rate W (Rayleigh's dissipation function; 
cf. Goldstein 1980) can be written 

W = F qi qi- (4-2) 
Therefore, to calculate T qi , we need to evaluate W and express it in terms of qi and qi. 

4.1 Viscous Dissipation 

The dissipation rate due to shear viscosity is given by (cf. Landau & Lifshitz 1987) 

W = - Ja ij u iJ d s x, (4.3) 
where ui is the fluid velocity and cry is the viscous stress tensor, 



cry = V [ iHj + - 3% V ' u J ' ( 4 - 4 ) 
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We denote by rj = pv the dynamical shear viscosity, where v is the kinematic shear viscosity. 
The bulk viscosity will be neglected. From equations (2.2)-(2.4), we have 



a\ a-i 
= — , u 2 ,2 = — , 

Mi. 2 = — A - S2, 

a 2 



1-3 



a 2 

= 0, otherwise. 



U 2 A = A + fi, 

ai 



Thus the viscous dissipation rate is given by 



o 3 

^3,3 = — , 
a 3 



(4.5) 
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03 
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(4.6) 



where v is the mass-averaged shear viscosity 



v = — I v dm. 

M 



(4.7) 



Since W in equation (4.6) is quadratic in dj, from equation (4.2), the dissipative forces 
are given by 

To* = 



Thus we have 



2 dqi 



(4.8) 
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(4.9) 
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oia 2 



The statement T$ = simply indicates that viscous forces conserve angular momentum 
(dJ s /dt = J-'ff, = 0), while they do not conserve fluid circulation since dC/dt = 7^ 0. 
In the presence of viscous dissipation, the dynamical equations (2.19)-(2.23) become 



01 



, , 10 _ /2di 
{•••}- - — v 



3k u 



Oi 



02 

«2 



ds\ 1_ 

03/ Oi 



(4.10) 
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{•• •} - - — v 

^(aiO-a 2 A) = {•••}- —p , 
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(-a 2 + aiA) = {•••} 



ai 
ai 

ai 
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a 3 / a 2 
a 2 y a 3 
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af — a, . 
-A, 

a\a 2 
_ af - a\ . 

v 2~^> 

aia^ 



(4.11) 
(4.12) 
(4.13) 
(4.14) 



where {• • •} denote the terms that already exist in equations (2.19)-(2.23) (This notation 
will be used throughout the paper). For the pressure term, expression (2.30) still applies 
here, but expression (2.34) must be modified as 




{•••} + io* £3 



a- 



(n = 0). 



(4.15) 



Equations (2.24)-(2.25) become 
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a\ 


CLi 


a 2 


a 2 




ai 


a 2 



, , 10 a, - a, 
{•••}+—* 



2 2 
ft-j^ ^ 2 



2 1 
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{•••} + — i>- " 

K n aiG2 



a 



: + ai 



A 



(4.16) 



Recall that the enegy dissipation rate is simply £ = W. For a quasi-static ellipsoid, this 
rate reduces to equation (6.4) in LRS4. 

4.2 Gravitational Radiation Reaction 

In the weak-field, slow-motion regime of general relativity, the emission of gravitational 
waves induces a back-reaction scalar potential $ re act which can be written as (Misner, 
Thorne & Wheeler 1970) 

G i(5) 
5 C 5 ij 



react 



(4.17) 



where the superscript (5) indicates the fifth time derivative and c is the speed of light. 
Here we choose X{ (i = 1, 2,3) to be the Cartesian coordinates of a fluid element in the 
instantaneous corotating frame (the body frame, with basis vectors along the principal 
axes). Then iff is the fifth time derivative of the reduced quadrupole moment tensor of 



(5) 

the system in the inertial frame projected onto the body frame. Expressions for are 
given in Appendix A. The radiation reaction force per unit mass on the fluid is — V<E> reac t. 
The energy dissipation rate is thus given by 



W = - / U • VQreact dm. 



(4.18) 
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Using equations (2.2)-(2.4), we get 

W = - Q««M^ |j [i { Sa iai + l { £a 2 a 2 + l^a 3 a 3 + lf 2 } 0(a? - a|)] , (4.19) 
where we have used 

JxiXjdm = Iij = -K n Maf5ij. (4.20) 

Since W here is linear in the dissipative forces due to gravitational radiation are 
given by (compare eq. [4.8]) 

From equations (4.19) and (4.21) we obtain 

(1 \ 2,G 
5 KnM ) 5^2? ° 2 ' 



5~ n / 5c* 33 



J> = o. 

Since JF^, = 0, we see that gravitational radiation reaction conserves the fluid circulation, 
i.e., dC/dt = = 0. But since T$ ^ 0, the total angular momentum is not conserved in 
general. Indeed, gravitational waves can carry away angular momentum as well as energy 
when the system is not axisymmetric. 

Including the gravitational radiation reaction, the dynamical equations (2.19)-(2.23) 
become 

a 1 = {...}-^Sa 1 , (4.23) 
a 2 = {•••}- |^4? a 2 , (4.24) 
a3 = {"-}-^4 5 3 ) a3, (4.25) 
|(a 1 0-a 2 A) = {---}-g4 5 2 ) a 1 , (4.26) 
j f (-a 2 Q + oiA) = {•••}- f^W (4.27) 

Since I^ +4^ + 4^ = 0, it can be shown easily that the expression for 2P C / p c is not 
affected by the presence of gravitational radiation reaction, i.e., equations (2.30) and (2.34) 
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still apply. Equations (2.24)-(2.25) are modified as 



a 2 






Cl2 




a\ 


Ol 


0-2 



1 \- ■ ■} + f- + - 

5c 5 \a,2 a\ 



4G T ( 5 ) 



(4.28) 



As before the energy loss rate is simply S = W. 



4.3 Secular Evolution of a Compressible Ellipsoid 

The secular evolution of incompressible ellipsoids in the presence of dissipation has 
been studied previously. Using a linear perturbation analysis, Roberts & Stewartson (1963) 
and Chandrasekhar (1970) demonstrated that the presence of either viscosity or gravita- 
tional radiation reaction induces a secular instability of incompressible spheroids on the 
Maclaurin sequence beyond the point where the Jacobi and Dedekind sequences branch 
off. Press & Teukolsky (1973) have studied the viscous evolution of a secularly unstable in- 
compressible Maclaurin spheroid, and Miller (1974) has integrated the Riemann-Lebovitz 
equations for incompressible ellipsoids including the effects of gravitational radiation reac- 
tion. 

Our dynamical equations represent a generalization to compressible ellipsoids of the 
equations considered by Press & Teukolsky (1973) and Miller (1974). We have integrated 
our equations for a variety of compressible ellipsoidal configurations, and have found results 
similar to those obtained in these previous studies for incompressible ellipsoids. 

The secular evolutionary paths of the ellipsoids can be understood clearly in terms of 
the energetics and conservation principles: dissipative forces always drive a system toward 
a state with lower energy, along quasi-equilibrium paths that hold conserved quantities 
fixed. Consider first the evolution of an ellipsoid under gravitational radiation reaction. 
Since the radiation reaction forces conserve the fluid circulation (cf. eq.[4.22]), such an 
evolution is always along a constant-C sequence. In Figure 1, we show the variation of 
the total energy along various constant-C sequences (with n = 1) as a function of the axis 
ratio a^/ai- The curves were obtained numerically by solving the Riemann-S equilibrium 
equations (LRS1, §5). Note that a given Maclaurin spheroid corresponds to a unique value 
of C, but two different constant-C sequences branch off: one sequence is Jacobi-like, with 
|0| > |A|, the other is Dedekind-like, with |A| > |0|. 4 For given values of C and a 2 /ai, 
the Jacobi-like configuration has higher energy than the Dedekind-like configuration. Also 
shown in Figure 1 is the Dedekind sequence (of configurations with A/O = oo) and the 
Jacobi sequence (A/O = 0), which bifurcate from the Maclaurin sequence ((12/ai = 1) at 



4 Such a characterization of the two branches (as used by Detweiler & Lindblom 1977) 
is not exact. When the deformation is very large, e.g., a 2 /ai < 0.1, we find that both 
configurations have |0| > |A|. However, except for such highly deformed configurations, 
the two branches do have |0| > |A| or |A| > |0|. More appropriate is using |C/0| < 2 to 
define Jacobi-like and |C/^| > 2 to define Dedekind-like, where ( is the vorticity in the 
corotating frame. 
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the point where the ratio T s /| W| = 0.1375 (in our approximation, this value is independent 
of n; see LRS1, §4). Note that the Dedekind sequence and the Jacobi sequence are adjoints 
of each other, and therefore have the same energy for given aija\ (see LRS1, §5.1). 

We see from Figure 1 that there exists a critical value of \C\ = \C sec \, equal to the abso- 
lute value of the circulation of a Maclaurin spheroid at the bifurcation point, below which 
a constant-C sequence has a Maclaurin spheroid as its minimum-energy state, and above 
which its minimum-energy state is a Dedekind ellipsoid. Therefore, an initial configuration 
with |C| < \C sec \ will evolve to become a Maclaurin spheroid under gravitational radiation 
reaction. For an initial configuration with |C| > \C sec \, the final state is a Dedekind ellip- 
soid. In this latter case, however, if \Q\ > |A|, the system will first evolve to a Maclaurin 
spheroid; only when some additional perturbation (e.g., viscosity) triggers the secularly 
unstable bar mode of the Maclaurin spheroid does the system evolve past the Maclaurin 
sequence, toward a final Dedekind ellipsoid. Note that the evolutionary timescale of the 
Dedekind-like phase is much longer than that of the Jacobi-like phase, as a result of the 
steep power-law dependence of the energy dissipation rate on O (cf. Lai & Shapiro 1994). 

The evolution of ellipsoids under pure viscous dissipation can be understood similarly. 
As viscosity conserves J (see eq. [4.9]) while dissipating energy, the evolution is along a 
constant-J sequence. From Dedekind's theorem (Ch69, Chap. 3; LRS1, §5.1), we know 
that the energy curves of constant— J sequences are identical to those of constant-C curves, 
e.g., a Jacobi-like constant-J sequence is adjoint to a Dedekind-like constant-C sequence 
with C = — J. Therefore, Figure 1 can be used again here, but after switching the curves for 
Jacobi-like sequences and Dedekind-like sequences. We see that the final state is either a 
secularly stable Maclaurin spheroid (when J < J sec ) or a Jacobi ellipsoid (when J > J sec )- 

The difference in the evolution for different polytropic indices n is illustrated in Fig- 
ure 2, where we show the E vs J curves of Maclaurin, Jacobi and Dedekind sequences 
for various n. In the presence of viscosity, an ellipsoid evolves vertically downward in this 
diagram, and the evolution terminates at a corresponding Maclaurin spheroid or Jacobi 
ellipsoid. We see that for small n <^ 1.5, a Dedekind ellipsoid evolves to a secularly stable 
Maclaurin spheroid; but for large n J> 1.5, a Dedekind ellipsoid first evolves toward a 
secularly unstable Maclaurin spheroid, and finally to a Jacobi ellipsoid. This qualitative 
difference is easy to understand: a highly compressible configuration can expand appre- 
ciably when deformed, giving a larger angular momentum (this can also be seen using the 
scaling relation for Riemann-S ellipsoids; cf. eq.[3.27] in LRS1). 

Similarly, Figure 2 can be used to understand the evolution driven by gravitational 
radiation reaction. In this case, the horizontal axis represents — C, and curves for Jacobi and 
Dedekind sequences are switched. Thus we see that for small n, a Jacobi ellipsoid evolves 
to a Maclaurin spheroid, while for large n, its final state is a Dedekind ellipsoid. One of the 
outstanding problems of 3D hydrodynamics with gravitational radiation is to verify this 
behavior using the exact hydrodynamic equations. The result has important consequences 
for the fate of nonaxisymmetric, rapidly rotating neutron stars, and for coalescing binary 
neutron stars (see Rasio & Shapiro 1994, §4.1). 

Lindblom & Detweiler (1977) have considered the combined effects of gravitational 
radiation reaction and viscosity on the stability of incompressible Maclaurin spheroids. 
Based on a linear analysis, they showed that when operating together, the two effects tend 
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to cancel each other. Evolutionary tracks of general ellipsoids for various ratios of the 
viscous timescale and the radiation-reaction timescale have been obtained by Detweiler 
& Lindblom (1977). Using our dynamical equations we can now extend this study to 
compressible ellipsoids (Lai & Shapiro 1994). 

4.4 Secular Instability Growth Time 

We have derived the secular instability growth time for a compressible Maclaurin 
spheroid, both with viscosity and gravitational radiation reaction (Lai & Shapiro 1994). 
The viscous instability growth time r v i s is given by 

while the growth time tqw of the instability driven by gravitational radiation reaction is 
given by 

t G w = — ■ ( 4 - 3 °) 

Here a is the eigenfrequency of the toroidal mode in the absence of dissipation, given by 
equation (3.30), and K n is defined by equation (2.7). In the the n = limit, these results 
agree with those given by Ch69 (§37) for the viscous instability and by Chandrasekhar 
(1970) for gravitational radiation reaction. 

A more complete discussion of the instabilities and secular evolution of rotating stars 
will be presented elsewhere (Lai & Shapiro 1994), together with an application to the 
emission of gravitational waves during core collapse. 

5. DYNAMICAL EQUATIONS FOR BINARIES: 
ROCHE-RIEMANN SYSTEMS 

Having studied in detail the evolution of single stars modeled as compressible Riemann- 
S ellipsoids, we now turn to binary systems. We will derive the general dynamical equations 
for a bound or unbound system containing a compressible Riemann-S ellipsoid and a point 
mass (a Roche- Riemann binary system). For parabolic orbits, Nduka (1971) first derived 
the dynamical Riemann-Lebovitz equations in the incompressible limit 5 . In addition to 
providing a generalization to compressible ellipsoids, our equations also determine the 
orbital dynamics self-consistently, thus allowing for general binary orbits. Equilibrium 
Roche-Riemann binaries in circular orbit have been studied in LRS1 (§8). 

In addition to the coordinates {a^, cf>, ifj} used in §2 to describe the structure of a single 
ellipsoid, we need to introduce new variables to specify the orbital motion and the relative 
orientation of the ellipsoid. For simplicity we consider only orbits in the equatorial plane of 
the ellipsoid (i.e., with the orbital angular momentum and the spin of the ellipsoid aligned 
along e 3 ). Two new coordinates are then needed to describe the orbit: a radial coordinate 



5 Misprints and errors in the original paper by Nduka have been pointed out by Luminet 
& Carter (1986; p. 224) and by Kosovichev & Novikov (1992). 
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r which measures the separation between the center of mass of the ellipsoid and the point 
mass, and an angular coordinate 9 which we take to be the true anomaly of the point mass. 
We also define a misalignment angle a = 9 — <fi between the axis a\ of the ellipsoid and 
the line joining the centers of the two bodies (see Figure 3) . The point mass is called M' 
and, following Ch69 and LRS, we denote the mass ratio by p = M/M' . 

The Lagrangian of the system is simply the sum of the single ellipsoid contribution, 
L s , given by equations (2.8)-(2.12), and an orbital contribution 7 r6, 

L = L s + L or b, (5.1) 

where 

Lorb = \nf 2 + ^/jr 2 9 2 - Wi. (5.2) 

The first two terms in equation (5.2) give the orbital kinetic energy, and the last term Wi 
is the interaction energy between M and M', given by 

W l = -GM' Jd 3 x ' ( 5 - 3 ) 

where p(x) is the density distribution within M. To quadrupole order, we have 

1 1 x f 1 



r 



2 



+ 2^[3(x.*) 2 -x 2 ], 



where f is the unit vector connecting M to M'. Since J(x • r)dm = 0, we have 

GMM' GM' T . 

Wi = ^3- (3/ rr - In - 7 22 - 7 33 ), (5-5) 

where 

I rr = d 3 xp(x)(-x ■ f) 2 = In cos 2 a + 7 2 2 sin 2 a, (5.6) 



and In is defined by equation (4.20). Combining (5.5) and (5.6), we obtain 
GMM' CM' 

Wi = ^[ /n (3cos 2 a - 1) + 7 22 (3sin 2 a -I)- 7 33 ]. (5.7) 

r Zr° 

Given the Lagrangian (eqs. [5.1] and [5.2]), the dynamical equations can then be 
obtained from the Euler-Lagrange equations (2.13). Of particular interest are the equations 
for <j> and 9. For q t = </>, we get 

dJ s - 3GM' ^ « /r nN 

=M= -^3-sin2a(/u -7 22 ), (5.8) 

where J s is the "spin" angular momentum of the star, given by equation (2.17), and Af is 
the tidal torque exerted on the star. For qi = 9, equation (2.13) gives 
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where J or b = pr 2 6 is the orbital angular momentum. Thus we see that the total angular 
momentum, 

J = J s + J orb , (5.10) 

is conserved, as expected. Equation (2.13) for qi = ip indicates that the fluid circulation 
given by equation (2.18) is also conserved, since tidal forces conserve the fluid vorticity. 

The complete dynamical equations for a Roche-Riemann system in the absence of 
dissipation can be written as 

d x = ai O 2 + A 2 - 2a 2 OA a 1 A 1 p + — - — 1 5— -(3 cos 2 a - 1), 

q n \nK n pcj ai r 6 

(5.11) 

a 2 = a 2 (Q 2 + A 2 ) - 2 fll OA - —a 2 A 2 p + ( - + ^^(3 sin 2 a - 1), 

q n \nK n pcj a 2 H 

(5.12) 

2tvG -^[5h P c \ 1 GM'a 3 
a 3 = a 3 A 3/ 9+ = — , (5.13) 

q n \nK n p c J a 3 r 6 

d , »n • • » 3GM'ai . n 

— (01O - a 2 A) = -aifi + a 2 A H — -= — sin2a, (5.14) 

at 2r 6 

d , . ./->.* 3GM'a 2 . . . 

— (— a 2 SZ + aiA) = a 2 S2 — aiA H = — sin 2a, (5.15) 

at 2r 6 

.. - 2 g(M + M') 3« n G (M + MQ r 2 2 1^2,0-2 n 21 

r = rO 5 — -7 a^(3cos a — 1) + a^(3 sin a — 1) — a 3 , 

r 10 r 

(5.16) 

•• 2r6 3 Kn G (M + MQ 2 2 , . 9 ... 
= — — ^ (o 1 -a 2 )sm2Q, (5.17) 

where # = O or ^, = O. For numerical integrations, equations (5.14)-(5.15) can be rewrit- 
ten as 
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(5.18) 
(5.19) 



The pressure term in equations (5.11)-(5.13) can be handled in the same way as for a 
single star. It is easy to verify that equations (2.30) and (2.34) still apply. 
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6. DYNAMICAL INSTABILITY OF ROCHE-RIEMANN BINARIES 

It is straightforward to show that, for an equilibrium system (i.e., when hi = = di 
and a = 0), the dynamical equations (5.11)-(5.19) reduce to the equilibrium equations 
for Roche- Riemann binaries derived and solved in LRS1 (§8.1). We now examine the 
dynamical stability of the equilibrium solutions. 

To study small dynamical oscillations, one could linearize the dynamical equations 
and calculate all the eigenfrequencies. This involves extensive algebra. However, we 
can show that the onset of dynamical instability is determined by a condition simi- 
lar to equation (3.8), with an appropriately constructed energy function. For this pur- 
pose, it is convenient to use a instead of 6 as an independent variable. Thus we define 
{li} = { a i? a 2, «3, r, a, (f>, ifj} here. The Lagrangian of the system as discussed in §5 can be 
written as 

i 1 • • 2 • • 

-K n M{h\ + hj + hi) + -I{4> 2 + V> 2 ) - -K n a!a2<j)ip 



10 

+ ^r 2 + ^r 2 (0 + a) 2 — U — W — Wi 



(6.1) 



The canonical momenta associated with the {qi} are 

P ai = -K n Mhi, P r = nr, P a = iir 2 (<p + a), P^ = J, P^=C. (6.2) 
o 

Because of the conservation of J and C, and also because a = for an equilibrium con- 
figuration, it is convenient to define a subset of the variables {cti} = {a±, a,2, 03, r}. The 
Hamiltonian can then be written as 

if (a,, a, P ai , P a , J, C) = 2(/ J^ /5) {Pi + Pi + Pi) + + «» P ^ J> C), (6.3) 

where 

E( ai ,a,P a ,J,C) = ^Pa + ^iJ-Pa+Cf + ^iJ-Pa-Cf + U + W + Wi, (6.4) 

(compare with eq. [3.3]), where I± is defined in equation (3.4). Hamilton's equations 
dP q Jdt = -dH/dqi for q t = a; give 

-KMd--(— ^ --(—} + a(—^\ (6 5) 

where for the second equality we have used dE/dP a = a, and we consider P a to be a 
function of ctj, J, C and a, and £ to be a function E(ai,a,a, J,C). Linearization of 
equation (6.5) about equilibrium, letting on = ai ;eq + Soti, yields 

5 Kn ttl ^ aj \daidaj) eq a \da l da) JC a \da i daJ JC ° V da i / 

52 ; 



Vjq f ^ a ( \ Ya{— 

^ <<y, \ da t da 3 J eq a \daida) J C a \da iy I C 
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(6.6) 



where the second equality follows because (d 2 E/ddida) oc sina (see eq. [5.8]) and \a\ <C 1 
near equilibrium. Similarly, Hamilton's equation for qi = r yields 

Now we substitute 5a; oc e %at in equations (6.6) and (6.7) and let a = 0. We obtain 

d 2 £ 



daidcti , 

4 J / eq 



0, for a = 0. (6.8) 



Clearly, in this expression, E should be evaluated at a = 0, the equilibrium value. Thus 
we can write E = E(oti; J,C) with fixed a = 0, and the onset of dynamical instability is 
determined from the condition 

d e M t; — q — ) =0, i, j/ = 1,2,... (onset of instability) (6.9) 



\daidajy eq 

where the partial derivatives are evaluated holding J, C fixed. This condition forms the 
basis of the stability analysis of binary systems presented in LRS1 and LRS4. Although 
the new degree of freedom associated with a was not introduced in our previous analyses, 
the explicit derivation given above shows that the stability condition is not affected. 

Alternatively, on quite general grounds, one can show that the stability conditions de- 
termined from equation (6.9) coincide with turning points along appropriately constructed 
equilibrium sequences (LRS1, §2.3). Specifically, a dynamical stability limit coincides with 
the point where the total equilibrium energy and angular momentum are both minimum 
along a sequence with constant circulation. Using both methods (eq. [6.9] and the turning 
point method), we have found in LRS1 (§9.2) that a Roche- Riemann binary can become 
dynamically unstable when the orbital separation is sufficiently small. This instability 
results from the strong tidal interaction, which can make the effective interaction potential 
between the two stars much steeper than 1/r, thereby destabilizing a circular orbit (cf. 
Goldstein 1980, §3-6; see also LRS2 for a qualitative discussion). The stability limits for 
various Roche-Riemann binary models have been tabulated in LRS1 (Tables 10 and 11). 

To illustrate how the instability develops, using our dynamical equations, we show 
in Figure 4 the time evolution of an unstable system with n = 1, p = M/M' = 1 and 
A = (corotating). The dynamical equations were integrated numerically using a standard 
fifth-order Runge-Kutta scheme with adaptive stepsize (Press et al 1992). At t = 0, an 
equilibrium solution is constructed for rja\ = 1.7, and r/R a = 2.406. This equilibrium 
solution is then perturbed by setting r = 10~ 3 (GM/Ro) 1 / 2 . For comparison, the results of 
an integration for a stable binary with r/a± = 1.8, r/R a = 2.427, and with the same applied 
perturbation is also shown. The dynamical stability limit along the corotating (Roche) 
sequence with n = 1 and p = 1 is at r/a± = 1.760, or r/R Q = 2.417. We see clearly in 
Figure 4 that, as the dynamical instability develops, a± increases while r decreases, and this 
is accompanied by the significant development of a tidal lag a > 0. Of course, the precise 
evolution of an unstable binary depends on how the initial configuration is perturbed. 
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7. TIDAL CAPTURE AND DISRUPTION OF A STAR 



BY A MASSIVE BODY 



Tidal interactions of stars and other fluid bodies have been discussed extensively in 
a number of different contexts. Fabian, Pringle and Rees (1975) originally proposed that 
tidal encounters between a neutron star and a main-sequence star might lead to the forma- 
tion of X-ray binaries in globular clusters (but see Rasio & Shapiro 1991; Kochanek 1992a; 
Rasio 1993). The first quantitative calculations of the orbital energy dissipation were per- 
formed by Press & Teukolsky (1977, hereafter PT) using a linear perturbation method. 
They were followed by many other studies using both linear theory (Lee & Ostriker 1986; 
McMillan, McDermott & Taam 1987; Kochanek 1992a and references therein) and numer- 
ical hydrodynamic calculations (Rasio & Shapiro 1991). The tidal disruption of stars by 
a massive black hole can provide a mechanism for fueling low-luminosity active galactic 
nuclei (AGN), and may also lead to observable flares in the luminosity of AGN. Many 
aspects of this problem have been considered previously (Hills 1975; Rees 1988; Evans & 
Kochanek 1989; Carter & Luminet 1985; Novikov, Pethick & Polnarev 1992), including 
relativistic effects (Laguna et al. 1993). Tidal disruption of small bodies (planetesimals) by 
protoplanets has also been discussed in the context of the Solar System formation (Boss, 
Cameron & Benz 1991; Sridhar & Tremaine 1992). 

The advantage of our ellipsoidal method for studying this problem is that it allows 
for the treatment of nonlinear effects, which are inevitable for close encounters. Using 
their affine stellar model, Carter and Luminet (1985, 1988) have studied in detail the 
tidal disruption of stars by massive black holes. Other studies using similar affine-type 
models include those of Kochanek (1992a), Kosovichev & Novikov (1992), and Sridhar & 
Tremaine (1992). There are several inconsistencies and differences in the results obtained 
in those previous studies. Thus we consider it worthwhile to re-examine the problem using 
our independent formulation of the dynamics, even though our equations are formally 
equivalent to those of the affine model. We include in some of our calculations the effects 
of fluid viscosity, which can be significant for the viscoelastic material in planetesimals 
(Sridhar & Tremaine 1992). 

For definiteness, we focus on the encounter of a fluid body of mass M with a point-like 
object of mass M' ^> M, and we consider only parabolic trajectories. A useful dimension- 
less parameter characterizing the encounter is 



where r p is the periastron separation. The quantity r\ is simply the ratio of the timescale 
for the periastron passage, ~ r p /v p , where v p is the velocity at the periastron, to the 

dynamical timescale of the star, ~ R^J 2 / (GM) 1 / 2 . When the relative velocity between 
M and M' at infinite separation is nonzero, the orbit is slightly hyperbolic. However, as 
long as v p ^> Woo, the trajectory remains very nearly parabolic close to periastron, where 
the tidal interaction is strongest. 




(7.1) 
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7.1 Dynamical Calculations for Compressible Ellipsoids 



For a given 77, we integrate equations (5.11)— (5.19) numerically to determine the dy- 
namical effects of the tidal interaction. Initially the star is placed on a parabolic orbit, 
far away from the massive body (r/R ^> ?7 2 / 3 /j _1 / 3 ), and it is assumed to be spherical 
(nonrotating). To avoid the singularity in equations (5.18)-(5.19) when a-y = 0,2, the initial 
configuration is slightly perturbed by increasing (decreasing) a\ (02) by ~ 0.1%. We have 
checked that the final numerical results (including the energy transfer) are independent of 
the exact amplitude of this initial perturbation. The total angular momentum, fluid circu- 
lation, and energy are conserved to very high accuracy throughout the evolution (typical 
error <> 10" 7 ). For n = 0, we also check the conservation of the volume: the product 
remains constant to within ~ 10~ 7 typically. For most calculations we use M'/M = 10 6 , 
but the results are essentially independent of the exact value as long as M'/M J> 10 3 . 

Typical results are illustrated in Figure 5 for an encounter with r\ = 2.8 and n = 1.5. 
Here t = corresponds to the time of periastron passage. After the encounter, the star 
becomes an oscillating Riemann-S ellipsoid with zero fluid circulation, but finite angular 
momentum. The angular momentum deposited in the star through the tidal interaction 
does not lead to uniform spin in the absence of fluid viscosity. 

7.2 Energy and Angular Momentum Transfer in the Linear Theory 

The energy transferred from the orbit to the star during a tidal encounter can be 
calculated exactly in the limit of linear perturbations using the theory developed by Press & 
Teukolsky (1977). Here we also calculate the corresponding transfer of angular momentum 
using the same formalism. 

The amount of energy transferred can be written as 

J3 „ _ d£ 



AE=- J dt j d 6 xp-±- VW, (7.2) 

where U is the gravitational potential of the point mass, U(x,t) = —GM'/\x. — r|, with r 
describing the orbital trajectory and x specifying the position of a fluid element inside M. 
The Lagrangian displacement £ of a fluid element in the star is assumed to be small, and 
can be decomposed into normal eigenmode components £ n / m , where {n, l,m} are indices 
specifying the eigenmodes. An important quantity in any discussion of tidal interactions 
is the dimensionless coefficient characterizing the coupling between the tidal potential and 
a particular normal mode (Zahn 1970; PT), 

Q nl = Jd 3 xpC nlm ■ V [r l Y lm {9A)\ = J*plr l+1 dr [C nl (r) + (I + l)^(r)] , (7.3) 

— * 

where, for spheroidal modes, we have expressed £, n im as a sum of radial and tangential 
components, 

imm(r) = [Cni(r)e r + r^ nl (r)V]Y lm (9, <P). (7.4) 
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The normal modes are normalized so that fd 3 xp£ n i m ■ i n >i'm' = 8nn'$ll'fimm'- The 
total energy transfer during a parabolic encounter can be written (PT) 

U ° 1=2 V T P ' 

where the dimensionless function Ti is given by 

I 

T l ( V )=2n 2 J2\Qm\ 2 Yl \ R nlm\ 2 - (7.6) 

n m=—l 

The function K n i m involves the eigenfrequencies u n i and is given in PT. 
The angular momentum transfer during an encounter is given by 

AJ Z = Jdt J d 3 x (p + Sp)r z = Jdt J d 3 x dp ("f^Q , ( 7 - 7 ) 

where dp = —V • (p£) is the Eulerian perturbation of the fluid density in the star, and 
t z = —dU/d(f) is the tidal torque per unit mass. Using a similar procedure as in PT, we 
obtain (Lai 1994b) 

AJ ^^r(^r) SU) ^ (7 - 8) 

where the dimensionless function Si is given by 

Stir,) = 2n 2 J2 l — E (-^)l^m| 2 . (7.9) 

n ra——l 

Note that K n i m is larger for m = —2, thus AJ Z in (7.8) is positive. Also note that 
contributions to AE and to A J from different terms in the sum are related by AJ n i m = 
{-m/u n i)AE n i m . 6 

In the incompressible limit (n = and T = Ti = oo), only f-modes of oscillation 
exist. These have eigenfrequencies ^/(GM/R 3 ) = 21(1 — I)/ (21 + 1) (see, e.g., Cox 
1980). For the dominant / = 2 quadrupole modes (including the Kelvin mode, cf. §3.3), 
ul 2 /(GM/R 3 ) = 4/5. The normalized eigenfunctions are £ r 02 = 2^ 02 = (8^/3)^^, for 
which we have Q 02 = (3 /2-k) 1 / 2 MR 2 . 



6 In eqs. (7.6) and (7.9), the index m does not correspond to the original mode index 
in Yira- Instead, contributions from m = 2 and m = — 2 modes have been re-grouped to 
derive these final expressions. In fact, it can be shown that the m = 2 and m = — 2 modes 
contribute to the tidal energy and angular momentum equally. 
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7.3 Large-?] Encounters: Comparison with Linear Theory 

In Figure 6, we compare our ellipsoidal results for the energy and angular momentum 
transfer (§7.1) with those of linear theory. In the limit where ?j > 1, the linear theory 
should be exact. The three values of n we have considered are n = 0, 1.5, and 2.5. For 
n = 0, one would expect the results from the two methods to agree precisely, since only 
the f-modes are involved in the linear theory, and the / = 2 quadrupole interaction is the 
leading order (the / = 2 mode structure obtained from the two methods is the same; see 
§3.3). We find that this is not the case, i.e., even in the limit of n — > oo, the results from 
the two methods do not agree. The same conclusion was reached by Kosovichev & Novikov 
(1992). Kochanek (1992a) considered only n = 1.5 and concluded that the discrepancy 
was probably due to small errors in the numerical integration of his affine-model equations. 
Kosovichev & Novikov (1992) attributed the discrepancy to slight deviations of the true 
f-mode displacement from an exact ellipsoid. They argued that in an ellipsoidal model, 
the stellar shape is too prolate, leading to stronger tidal interaction. We question the 
validity of this interpretation, since for a circular binary in equilibrium, it can be shown 
explicitly that the energy shifts obtained from the linear theory and from the ellipsoidal 
(Roche-Riemann) model are identical (Lai 1994b). Note, however, from Figure 6, that the 
difference in the amount of energy dissipated is less than 40%. The resulting tidal capture 
radius is therefore hardly affected, since the dependence of AE on r p is extremely steep. 

In the compressible cases (n > 0), other types of nonradial oscillation modes, which are 
absent in the ellipsoidal model, can be excited by the tidal interaction. As discussed in §2.1, 
the ellipsoidal model incorporates only the / = 2 f-modes of oscillation. In general, p-modes 
are also excited, and, if Ti > T, one should also include g-modes. However, the coupling 
coefficients Q n i (eq. [7.3]) of the p-modes and g-modes are typically much smaller than that 
of the f-modes (Lee & Ostriker 1986). When g-modes exist (Fi > T), their contribution 
to the tidal energy transfer is larger than that of the f-modes at large separation (see the 
dotted lines in Fig. 6). The reason is that at larger separation, "resonances" occur since 
the (very low) frequency of the tidal driving force is always close to that of some high-order 
g-mode. These resonances lead to a more efficient energy transfer at large separation. 

7.4 Close Encounters: The Tidal Disruption Limit 

As i] decreases, the amplitude of the tidal perturbations increases and linear theory 
eventually breaks down. Below the tidal disruption limit at some critical rj = rjdis, the 
amount of energy transferred to the star exceeds its original binding energy (the stellar 
energy E s , eq. [2.26], becomes positive), and the star is left unbound after the encounter. 
In Figure 7 we show the results of two calculations for n = 1.5. For one encounter we 
used 7] = 1.85, slightly larger than the disruption limit rjdis = 1-84, and for the other 
we used rj = 1.83 < rjdis- For rj > rjdis, the star remains bound after the periastron 
passage, although the axes experience large-amplitude nonlinear oscillations when rj is 
very close to rjdis- For rj < rjdis, the star becomes unbound (E s > 0); at least one of the 
axes keeps increasing monotonically in time after the encounter, leading to a progressively 
more and more elongated structure. This behavior is identical to that found by Kosovichev 
k Novikov (1992) for n = 0. 
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In Table 2, we list the values of rjdis for different polytropic indices. For n = 0, 
our results agree precisely with those of Kosovichev & Novikov (1992) and Sridhar & 
Tremaine (1992), but not with those of Luminet & Carter (1986). For larger n, r^ s 
decreases, since the tidal interaction is less effective for more centrally concentrated objects. 
For comparison, we also list in Table 2 the absolute Roche- Riemann limit, r) rr , which 
corresponds to the minimum separation for an equilibrium configuration, among all Roche- 
Riemann binaries, to exist in a circular orbit (see LRS1, §8.2). Note that the Roche limit 
for corotating binaries, as well as the irrotational Roche-Riemann limit for C = binaries 
(LRS3), all correspond to larger separations than rj rr . From Table 2 we see that for all 
cases, rj rr > Vdis, which is intuitively expected. 

We note that even when rj > r] dis: the object can become highly elongated (see Fig- 
ure 7), with ai ^> a2,a3, even though on simple energetic grounds it is still bound. In 
reality, such needle-like object is likely to be subject to the "sausage" instability of infinite 
cylinders (Chandrasekhar 1961), and would break up into small pieces. Of course, our 
ellipsoidal model is not capable of treating this process. 

7.5 Effects of Viscosity 

Viscous dissipation forces can be easily incorporated into the dynamical equations 
(5.11)-(5.19). Since the motion of the center of mass of the star is not affected by viscous 
dissipation (which depends only on the shear stresses inside the star), the viscous forces 
derived for an isolated star (§4.1) can be directly applied to binaries. The dynamical 
equations (5.11)-(5.15) are modified in exactly the same way as in equations (4.10)-(4.16) 
for single stars. Also, since T r = Tq = 0, equations (5.16)-(5.17) remain unchanged. We 
assume that V is a constant during the dynamical evolution. 

Typical results are illustrated in Figure 8. Here the fluid viscosity has the value 
v = 0.01 (GMRo) 1 / 2 . We find that the disruption limit rjdis is smaller than in the inviscid 
case. For rj > rjdis — 1.72, the star is still bound after the encounter (E 8 < 0). Comparing 
with Figure 7, we see that the fluid viscosity damps out the large-amplitude, nonlinear 
oscillations of the ellipsoid after the encounter. Since the fluid circulation is not conserved 
in the presence of viscosity, the fluid does not remain irrotational. Instead, because viscous 
forces tend to damp out the differential rotation (i.e., A — > after the encounter), the star 
evolves to become an equilibrium Jacobi ellipsoid on the viscous dissipation timescale. 
This is very different from the undamped, inviscid case. When the post-encounter stellar 
energy E s is positive, as in the case where rj = 1.71 in Figure 8, the star becomes unbound, 
with one of the axes increasing monotonically. This is similar to what we found in the 
inviscid case. 

In Figure 9, we show the disruption limits as a function of the fluid viscosity for three 
values of n. The maximum physical value for the viscosity is v (GMRo) 1 / 2 , correspond- 
ing to momentum transport across the whole star in a dynamical timescale. In all cases, 
the disruption limit of a viscous body is smaller than that of inviscid body, as a result of 
the viscous damping of kinetic energy. This result is in agreement with the conclusions 
reached by Sridhar & Tremaine (1992). Note, however, that the reduction of rjdis is signifi- 
cant only for very large viscosities, V J> 10~ 2 (GMRo) 1 / 2 . Typical viscosities in stars (even 
the large turbulent viscosity in convective envelopes) always remain <> 10~ 4 (GMRo) 1 / 2 . 
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8. SECULAR EVOLUTION OF BINARIES DRIVEN BY VISCOSITY 



As discussed in §7.5, the viscous dissipation forces can be easily incorporated into the 
dynamical equations for the binaries. In this section, we consider the secular evolution 
driven by viscous dissipation of a binary in a circular orbit. This subject has been well 
studied in the literature (see, e.g., Goldreich & Peale 1968 for a review). The basic ingre- 
dients were already established in the the weak friction theory developed by G. Darwin 
more than a century ago. Our purpose in this section is to consider some aspects of this 
classical astronomical problem in the context of our compressible ellipsoid model, which 
naturally extends the early treatments to the nonlinear regime. 

8.1 Tidal Lag Angle Due to Viscosity 

Our present understanding of the tidal evolution of binary systems is largely based on 
the weak friction model. In this model, the viscosity of the star induces a lag angle between 
the static tidal bulge and the direction to the companion. This results in a torque and 
angular momentum transfer between the spin of the star and the orbit. This mechanism 
can lead to synchronization of the spin with the orbital motion, and a corresponding 
evolution of the binary orbit. Any degree of nonsynchronization is necessarily associated 
with a tidal lag angle, given by 

An vRAn , 

where AO = O — fi s , uj ~ (GM/R 3 ) 1 / 2 is the fundamental frequency of the star, and 
tvisc ~ R 2 /v is the viscous dissipation time. 

We can easily derive the exact result in the limit of large binary separation using our 
ellipsoidal model. For large r/R we have a\ ~ a 2 so that J s ~ — C (see eqs. [2.17]-[2.18]). 
Using equation (5.8) for dJ s /dt and dC/dt = -uMK{a\ - a\) 2 / '(aia 2 ) 2 (see eq. [4.9]), we 
find 

. n . 10r 3 (a 2 -a 2 2 \ 

™ 2a ~ ^A Q 7TTjnr9 ( -L- 2 ) . (8.2) 

For large r we have 

o 

a\a 2 2 M \r 
(see eq. [A25] in LRS4), and the tidal lag becomes 

10 . i?AR 
^T (5 - n) GM- (8 ' 4) 

providing the coefficient of proportionality in equation (8.1) (recall that A ~ O — Q s ). 

8.2 Viscous Evolution 

In the presence of viscosity, only synchronized binaries can be in true equilibrium. 
Any degree of nonsynchronization necessarily leads to evolution. However, when the vis- 
cosity is sufficiently small, the binaries evolve slowly along a sequence of quasi-equilibrium 



3K n GM'R 2 V. a ± a 2 



a\-a\ 15 M' ( R\ 

^^nirri- , (8-3) 
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configurations. Thus, in most cases, dynamical calculations are not needed to follow the 
viscous evolution. Since the viscous forces conserve angular momentum, the binary evo- 
lution driven by viscosity proceeds along sequences of constant total angular momentum 
J. This has been discussed extensively in LRS4 (§6). Depending on the value of J, the 
final fate of the binary can be qualitatively different. Consider the example in Figure 10, 
where we show the variation of the total energy along three different equilibrium sequences: 
one is the corotating sequence, and the other two have constant J. All three sequences 
are for n = and M/M' = 1. A critical angular momentum J = J cr u is the minimum 
angular momentum along the corotating sequence. When J > J cr it, the constant- J curve 
intersects the corotating curve, at two points corresponding to the minimum or maximum 
of the energy along of the constant- J sequence. For J < J cr u, no configuration along the 
constant- J sequence lies on the corotating sequence. Clearly, for J > J cr u, the viscous 
forces drive the binary toward a stable corotating configuration, while for J < J cr it, the 
viscous forces drive orbital decay to final coalescence. More detailed discussions of these 
points can be found in LRS4. Values of J crit for general n and p can be found in Table 10 
of LRS1. 

When a dynamical instability is encountered during secular orbital decay, dynamical 
calculations are needed to follow the subsequent evolution. In Figure 11 we show an 
example of such a dynamical calculation. Here p = M/M' = 1, n = 1, and the viscosity 
v = 0.01 (GMRo) 1 / 2 is assumed to be constant throughout the evolution. The initial state 
is constructed by solving the equations for an equilibrium Roche-Riemann configuration 
with fn = — A(of + a|)/(naia2) = —4 and r/a\ =3 (LRS1, §8), corresponding to a total 
angular momentum J = 1.220. We set a = at t = 0. After a transient oscillation, 
a attains a small but finite value, corresponding to the tidal lag derived in §8.1. The 
reason is that, as mentioned before, in the presence of viscosity, the binary is not exactly 
in equilibrium state, and and viscosity-induced tidal lag is inevitable. During the secular 
evolution, the binary stays very close to an equilibrium evolution track with constant 
J. However, as the dynamical stability limit is approached, the evolution becomes much 
faster, and a large dynamical tidal lag develops. This dynamical behavior is similar to that 
shown in Figure 4. 

9. BINARY EVOLUTION DRIVEN BY GRAVITATIONAL 
RADIATION REACTION 

In this section, we incorporate gravitational radiation reaction in our dynamical equa- 
tions for binaries. This allows us to study binary coalescence driven by the emission of 
gravitational waves. In §9.2 we consider the coalescence of a Roche-Riemann binary. In 
particular, we study how the dynamical instability is approached as the binary orbit de- 
cays. Applications to coalescing neutron star binaries will be presented in a forthcoming 
paper (Lai, Rasio, & Shapiro 1994c). 

9.1 Dynamical Equations Including Gravitational Radiation Reaction 

To calculate the gravitational radiation reaction forces, we consider two coordinate 
systems (see Fig. 3): the body coordinate system, centered on M, with basis vectors {ei} 
spanning the principal axes, and the orbital coordinate system, centered at the CM of the 
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system, with basis vectors {e-}; ej is along the line joining M and M', e 2 is perpendicular 
to ej in the orbital plane, and e 3 is perpendicular to the orbital plane. The two coordinate 
systems are related by 

x\ = x\ cos a + x<i sin a — r cm , 

X2 = — x\ sin a + x 2 cos a, (9-1) 
£3 = x 3 . 

We now write the gravitational radiation back-reaction potential as 



'react = XjXj. (9-2) 



This has the same from as equation (4.17), except that here we have used a different 

(5) 

coordinate system: i~. is the fifth derivative of the reduced quadrupole moment tensor of 

the system projected onto the orbital frame. Expressions for l~. are derived in Appendix 
A. 13 
For M', the velocity is simply 

v = C e I + ftorbr' cm e2, (9.3) 

where r' cm is the distance between CM and M' . Thus the contribution to the dissipation 
function from M' is 

%' = -^[M'r' cm f' c J§ + M'r^JWg]. (9.4) 
The fluid velocity in M can be written as 

v = u + u orb = w^e; + (-f cm ej - fi orb r cm e 2 ), (9.5) 

where u is the fluid velocity relative to the CM of M and u or & is the orbital velocity. From 
equations (2.2)-(2.4), we have 

u-i = —xi + ( —A - O I x 2 , 



Therefore, 



ai \a 2 

„ 2 = ^x 2 +f-^A fO|.r, (9.0) 
a 2 V a i 

d 3 

u 3 = — x 3 . 
a 3 



1G 

V • V^react = [lifcCfc + U orb ] ■ — iffx^ej. (9.7) 



5c 5 J 1 

Using equation (9.5) and the relation between {e^} and {e-}, we get 
2G 



V • V^react = ^ 



u\x~j cos a — 1^ sin aj + U2XJ (l— 1 sin a + cos ctj 

i _i(5) _ • _t(5) _ q _t(5) 

-+- u 3 Xji^j r cm Xji^j r cm i i oro x ji^j 



(9.8) 
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Substituting equations (9.1) and (9.6) in this expression, integrating over the mass distribu- 
tion in M, and adding the contribution from M' (eq. [9.4]), we obtain the total dissipation 
rate 

W =W M + WM' = - / V • V^react dm + W M > 



2G fl 

5c^ 



-K n M 



fi(5) 2 , t(5) • 2 x(5) ■ \ 

I ijj cos a + sin a — ij 5 smza I aiai 
+ ^Ij^ sin 2 a + ig^ cos 2 a + ig^ sin 2a j 0202 + ^0,36,3 
+ ^lg } cos 2a + (ig } - lg } ) \ sin 2a^) (a 2 - a 2 2 )Q 



where we have used equation (4.20) and Mr 2 m + M'r' 2 m = [ir 2 . 
Using equation (4.21), the dissipative forces are then given by 



<Fa\ — 



2G fl 



3~ao — 



5c 5 

2G fl 



a 2 



5c 5 
2G 
5c* 



(^KnM^j (ijj cos 2 a + ig" 1 sin 2 a — ig* sin 2aj a\, 
K n Mj (i^ sin 2 a + ig^ cos 2 a + ig^ sin 2a j a 2 , 



2G fl 



5c 5 V5 



> n M 1 I ig } cos 2a + (i^ - 1^) ^ sin 2a ) {a{ - a%) 



(5) T (5), 



= 0, 



2G (5) 

(5) 2 



(9.9) 



(9.10) 



Therefore, the dynamical equations for binaries, equations (5.11)-(5.17) are modified as 



«!={•••} 

a 2 = {---} 



2G 
2G 

|(a 1 0-a 2 A) = {.--}-g 
_(- 02 fi + ai A) = {...}-^ 



lg } cos 2 a + lg } sin 2 a - ig sin 2a 

i(5) • 2 , t(5) 2 , t(5) • r, 

Ijj sin a + cos a + lj 2 sin 2a 



ai, 

a 2 , 



[g } cos 2a + - lg } ) sin 2a 



ig cos 2a + -(1^ - Ig ; ) sin 2a 
2 



r(5) x(5)x 



ai, 
a 2 , 



(9.11) 
(9.12) 
(9.13) 
(9.14) 

(9.15) 
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(9.16) 
(9.17) 



Since lg } + lg } + lg } = 0, the expression for 2P C / p c is not affected by the presence 
of gravitational radiation reaction, i.e., equations (2.30) and (2.34) still apply. Finally, 
equations (5.18) and (5.19) become 



= 


( a 2 


ai 








\ai 


a 2 


A = 


( a 2 


ai 




[at 


a 2 



-1 r 



{•••} + 
{•••} + 



2G 

5c^ 

AG 

5c^ 



lg ) cos2a+-(I^ I 



(5) 



I~ ) sin 2a 



a\ a 2 
— + — 

a 2 ai 



l(5) 

lj2 COS 



2a + -lg ) )sin2a 



(9.18) 

Note that the equation = again guarantees that the fluid circulation defined by 
equation (2.13) is conserved. 



9.2 Orbital Decay of Roche-Riemann Binaries 

With the dynamical equations derived in §9.1 and the expressions for l|p derived in 
Appendix A, we can now calculate the orbital evolution of a general Roche-Riemann binary 
driven by gravitational radiation reaction. In particular, we can study the dynamical 
behavior of a coalescing neutron-star-black-hole pair prior to final merging, at least when 
general relativistic effects are not too important (i.e., when tgr ~ 6G(M + M')/c 2 < 
(1 + M'/M) 1 / 3 ^; see LRS3, §5). 

At large separation, the orbital decay is secular, and the binary evolves along a Roche- 
Riemann equilibrium sequence with constant C (since gravitational radiation reaction forces 
conserve C). This quasi-static evolution has been discussed extensively in LRS3. At smaller 
r, when the dynamical stability limit is approached, the dynamical equations must be used. 
In Figure 12 we show an example of such a dynamical calculation for p = 1, n = 1, and 
R a = 5GM/ c 2 (a typical value for neutron stars; see LRS3) . The fluid viscosity is assumed 
to be zero. At t = 0, we construct an equilibrium Roche-Riemann configuration C = 
and r/ai = 5. Initially, the binary closely follows the equilibrium constant— C sequence. 
As the dynamical instability develops, both the radial velocity and the dynamical tidal 
lag increase considerably. Thereafter the two stars merge hydro dynamically in just a few 
orbits. This qualitative behavior has already been observed in the simplified calculations 
we presented in LRS2 and LRS3. The development of a large lag angle is not surprising. It 
arises because of the finite time necessary for the star to adjust its structure to the rapidly 
changing tidal potential (cf. Lai 1994a). 

The small but finite value of a observed in Figure 12 in the limit of large orbital 
separation can be calculated as follows. Including gravitational radiation reaction, the 
rate of change of the spin J s is given by 

^ = ^sin2a(/ 11 - I 22 ) (9.19) 
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(see eq. [5.8]), where T$ is given by equation (9.10). Using the expressions derived in Ap- 
pendix A, we see that the dominant contribution to T§ comes from the terms proportional 

to lg } ~ \m\ rh [ir 2 . Thus 



dJ, 3GM' 



dt ~ 2r 3 



sin2«(/H - 7 22 ) - |j Q««M^ (lg } cos 2a) (a 2 - a 2 2 ). (9.20) 



As in §8.2, consider the limit of large r, so that J s ~ —C. Since dC/dt = (cf. eq. [9.10]), 
we have dJ s /dt ~ 0. Equation (9.20) then gives 

64 fi^G^r 2 64 /GM\ 5/2 / M'\ 3/2 

Note that this result is independent of the equation of state (the polytropic index n does 
not appear). In Figure 13, we show the numerical results for a as a function of r. We see 
that the value of a is indeed given by equation (9.21) when r/R a ^> 1. However, at smaller 
separation, for r/R a <> 4, the lag angle can become considerably larger (by as much as an 
order of magnitude) than predicted by expression (9.21). This is a result of the dynamical 
instability. 
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APPENDIX A: Evaluation of li? 

13 

Consider a triaxial body with its diagonal moment of inertia tensor as defined in the 
body frame given by Uj (eq. [4.20]). We now calculate li^P (cf. eq. [9.2]), the components 

of 1^ projected onto a "projection frame." To define this frame, consider Figure 3, but 
for now neglect M' by focusing on M: the body coordinates are Xi and the projection 
coordinates are xj, but for now center the projection frame on O. In the special case of a 
single star, when the projection frame coincides with the body frame, the procedure for 
calculating i~ has been provided by Miller (1974). Below we generalize this procedure to 
the case where {x^ and {xj} are different, as we have for binaries (see §9). 

The coordinates {xi} and {xj} are related to the inertial coordinates by {Xi} by 

£, = T ia (</>)X Q , x-i = T- ia (0)X a , (Al) 

where 

T(</>) = I -sin^ cos<£ I , (A2) 




and similarly for T(6>). The components of the reduced quadrupole tensor in the inertial 
frame are 



Thus we have 



d 5 

/ d 5-m \ d 



= T. a (e)T ]p (e)^[Ti (0)4 (0 )Ifc/ 



5 / j5-m \ m 

_ \ ' ^5 / " 1 \ \ "~ (~im r>P T>m-P 

~ 2^ C m [jtt5=^ ikl J 2^ C P H ik H jl ' 



where the constants are binomial coefficients, and 



In general, equations (A4)-(A5) are complicated to evaluate. In the quasi-static limit, 
when \dcii/dt\ <sC \Clai\, simple expressions for if) can be derived. To lowest order in dcii/dt, 
we have 

f -E J **E c l R l R ]7 + 5 E ^ E m 

k p=0 k p=0 

where we have used 

tf = tImtUWm - \hk5 a ^ {Al) 
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and dlkk/dt ~ in the quasi-static approximation. 

The matrix R p can be now evaluated, keeping terms up to order Cl. After some 
algebra, we obtain 



ij 



sin 2a cos 2a 
=160 5 (/n - I 22 ) | cos2a -sin2a 















cos 2a 



sin 2a 



+ 40O 3 [O(/n -I 22 ) +2fl(I n -I 22 )} -sin2a -cos2a 















(A8) 



(5) 

We can now write down the components that appear in our dynamical equations. 
For the single star case, the projection frame coincides with the body frame. Setting a = 
in equation (A8), we see that the only nonzero components are 



l(5) 
1 11 _ 



r(5) _ 



L 22 



= 40O J [O(/ n - I 22 ) + 2Q(I n - 7 22 )], 



l<g =1$ = 16tf(I n -I 22 ), 



(A9) 



where the la are defined in equation (4.20). 

To evaluate for a Roche- Riemann binary, we now move the origin of the projection 
coordinates {xj} back to the system CM as in Figure 3 and equation (9.1). But then we 
can decompose the total moment of inertia into two contributions: the orbital part due to 
two point masses M and M', and the fluid part due to the ellipsoid M. The contribution 
from the ellipsoid is given by equation (A8). The resulting nonzero components of are 

l ii = -4? = 16ft 5 (In - I22) sin 2a + Wn 3 orb [2n orb /jrf + 2Q or6/ ur 2 ] 

+ 40O 3 [O(j n - / 22 ) + 2f2(/n - I 22 )] cos2a, 

lg } = lg } = 16fto r &^ 2 + 160 5 (/n - I22) cos 2a 

- 40n 3 [n(in - J 22 ) + 2f2(/n - I 22 )\ sin 2a. 



These expressions obviously do not apply to certain dynamical situations, such as 
the hyperbolic fly-by of a star past a black-hole. In such a case the variables r and 
change rapidly and higher-order time derivatives must be retained. For two point masses, 
the expressions for I~ J can be evaluated analytically with repeated use of the equation of 
motion. 
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TABLE 1 



Radial Pulsations of Polytropes" 



n 


O 


<7r, 1 


0. 


1. 


1. 


0.5 


1.355 


1.364 


1.0 


1.877 


1.913 


1.5 


2.706 


2.793 


2.0 


4.137 


4.305 


2.5 


6.909 


7.155 


3.0 


13.27 


13.27 



NOTE: a Here = ^/[(3r - 4)GM/R 3 ] and T = Ti = 1 + 1/n. 

TABLE 2 



Disruption Limits and Roche- Riemann Limits' 1 



n 


Vdis 


Tj rr 


0. 


2.21 


3.806 


0.5 


2.08 


3.427 


1.0 


1.96 


3.054 


1.5 


1.84 


2.698 


2.0 


1.71 


2.368 


2.5 


1.584 


2.117 



NOTE: a Here rjdis is the disruption limit for encounters between a star and a 
massive body, r] rr is the Roche-Riemann limit for equilibrium binaries in circular 
orbits. 
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FIGURE CAPTIONS 



FIG. 1. — Secular evolution tracks of a Riemann-S ellipsoid with n — 1 driven by gravita- 
tional radiation reaction. The energy of an ellipsoid as a function of the axis ratio aija\ is 
shown along various equilibrium sequences. Dedekind-like sequences (|A| > are shown 
on the left, Jacobi-like > |A|) on the right. The thick solid line corresponds to the 
Dedekind and Jacobi sequences, while the other lines correspond to constant-C sequences: 
C = C/(GM 3 J R ) 1 _/ 2 = -0.48 (dotted line), C = -0.4 (dashed line), C = -0.32 (long 
dashed line) and C = —0.25 (dotted-dashed line). The solid round dot marks the point 
of bifurcation. This figure can also be applied to pure viscous evolution, in which case 
the curves for Jacobi-like sequences and Dedekind-like sequences are switched, and C is 
replaced by J = —C =constant. 

FIG. 2. — Energy as a function of angular momentum for Maclaurin (solid lines), Jacobi 
(dotted ines) and Dedekind (dashed lines) equilibrium sequences, for n = 0, 1, 2. 

FIG. 3. — A sketch of the coordinate systems used for binaries. 

FIG. 4. — Evolution of the binary separation r, the axes a^, and the lag angle a for 
Roche-Riemann binaries, with p = M/M' = 1 and n = 1. The initial configurations are in 
equilibrium, and corotating. The left panels show the evolution of a dynamically unstable 
binary (r/a\ = 1.7, r/R a = 2.406), while the right panels show that of a stable binary 
(r/ai = 1.8, r/R a = 2.427). Both equilibrium configurations are perturbed by imposing a 
small radial velocity. 

FIG. 5. — Evolution of the axes, energy, and angular momentum of a star with n = 1.5 
during a parabolic encounter with a massive (M' ^> M) body. The parameter rj = 2.8 (eq. 
[7.1]). The solid line in the middle shows the total (conserved) energy in the system. 

FIG. 6. — Energy transfer and angular momentum transfer to the star during its parabolic 
encounter with a massive body. The heavy lines are the results obtained by integrating 
numerically our equations for a compressible ellipsoid, while the lighter lines show the 
results of linear perturbation theory (eqs. [7.5] and [7.8]). The solid lines are for n = 0, 
the short-dashed lines for n = 1.5, and the long-dashed lines for n = 2.5. In the linear 
theory, all modes are calculated assuming Ti = T, except in the case of n = 2.5, where we 
have also included curves (dotted lines) showing the results when Ti = 5/3. 

FIG. 7. — Same as Fig. 5, but for r\ = 1.85 (left) and rj = 1.83 (right). The dashed lines 
in the middle show the orbital energy. 

FIG. 8. — Same as Fig. 7, but here the viscosity is nonzero, with v = 0.01(GMR o ) 1 / 2 . 
On the left r\ = 1.73 and on the right r\ = 1.71. 

FIG. 9. — The disruption limit r/dis for the parabolic encounter of a star with a massive 
body, as a function of the fluid viscosity v. The viscosity is assumed to be constant during 
the encounter. The three curves correspond to n = 0.1 (solid line), n = 1.5 (short-dashed 
line) and n = 2.5 (long-dashed line). 

FIG. 10. — The total energy E eq of several equilibrium Roche-Riemann sequences for 
p = M/M' = 1 and n = 0, as a function of the binary separation r. The solid line is 
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for the corotating sequence (A = 0), the dotted line for a sequence with constant J = 
J/(GM 3 R ) 1 / 2 = 1.4 > J crit = 1.375, and the dashed line for a sequence with constant 

J = 1-3 < Jcrit- 

FIG. 11. — The evolution of a Roche-Riemann binary driven by viscosity. Here p = 1, 
n = 1, and z/ = 0.01 (GMRo) 1 / 2 . The initial configuration is an equilibrium Roche- 
Riemann binary with fn = —4 and r/ai = 3, corresponding to J = l.220(GM 3 R o ) 1 / 2 . 
Here r is the binary separation, ai are the ellipsoid's axes, A measures the degree of 
nonsynchronization (A ~ O — fi s ), a is the tidal lag angle, and E is the total energy of the 
system. 

FIG. 12. — The evolution of a Roche-Riemann binary driven by gravitational radiation. 
Here p=l, n=l, v = and R a = 5GM/c 2 . The initial configuration is an equilibrium 
Roche-Riemann binary with fn = —2 and r/a-y = 5. All quantities are defined as in Fig. 11, 
and v r = f is the radial velocity. The axes for an equilibrium irrotational Roche-Riemann 
sequence are shown as lighter lines, while the dynamical values are shown as thicker lines. 
The total energy for an equilibrium sequence is also shown as a dotted line. 

FIG. 13. — The tidal lag angle as a function of binary separation r. Here p = 1, n = 0.5 
and R a = 5GM/c 2 . The solid line is from our dynamical calculation, the dotted line is an 
analytic expression for large r (eq. [9.21]). 
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